
#######
###define function to estimate 6 types of models at various bandwidths
#######

library(lmtest)
library(sandwich)
library(date)
library(lubridate)
library(doBy)




models<-function(data, N, dv,memo=NULL, placebo=F, model=NULL, int=NULL, int.var=NULL, cluster=NULL, clust.var=NULL){
	
	cut2<-memo
	
	
	if(model%in%c("ols",NULL) & int!="yes"){
	
	if(placebo==T){
	cut2<-memo
	
data$period<-NA
data$period[data$dayno<cut2]<-"pre-memo"
data$period[data$dayno>=cut2 ]<-"post-memo"
data$period<-factor(data$period, levels=c("pre-memo", "post-memo"))
table(data$period)
	
	data$distance<-data$dayno-cut2
	table(data$distance)

	}
	
form<-as.formula(paste(dv,"~","dayno"))

means<-summaryBy(form, data=data, FUN=mean, na.rm=T)##get daily hit rates
means$lag.hit<-NA
means<-means[order(means$dayno, decreasing=F ),]




for(i in 2:length(means$dayno)){
	
	means$lag.hit[i]<-means[ i-1  ,   2]
	
}

names(means)[3]<-"lag.hit.new"


data<-merge(data, means, by="dayno")
head(data)
		
	
	data<-data[order(data$distance), ]
model.names<-c("loc. linear","loc. linear squared","loc. linear + lagged dv", "squared + lagged dv","cubic","diff")
results<-list()

coefs<-NA
ses<-NA
se.std<-NA


for(i in min(N):length(N)){
	
	results[[i]]<-cbind.data.frame(model.name=model.names, coef=coefs, se=ses, se.std=se.std)
}


names(results)<-N

if(cluster==T){


for(j in min(N):length(N)){
	
	if(N[j]>=5){
	
temp<-na.omit(data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),c("period","distance","weapon.hit","lag.hit.new","dayno",clust.var)])
dim(temp)	

	temp<-temp[order(temp$distance), ]

	##local linear
m<-ols(temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),dv]~period+distance+period*distance, data=temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),],model=TRUE, x=TRUE, y=TRUE)

m2<-lm(temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),dv]~period+distance+period*distance, data=temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),])


if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear"]<-NA
results[[j]]$se[results[[j]]$model.name=="loc. linear"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="loc. linear"]<-NA
	}

if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear"]<-m$coefficients["period=post-memo"]
new<-robcov(m, cluster=temp[,clust.var])
results[[j]]$se[results[[j]]$model.name=="loc. linear"]<-sqrt(vcov(new)["period=post-memo","period=post-memo"])
results[[j]]$se.std[results[[j]]$model.name=="loc. linear"]<-summary(m2)$coefficients["periodpost-memo","Std. Error"]
}

##local linear squared
temp$distance2<-I(temp$distance^2)
temp$distance3<-I(temp$distance^3)

m<-ols(temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),dv]~period+distance+period*distance + distance2+period*distance2, data=temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),],model=TRUE, x=TRUE, y=TRUE)

m2<-lm(temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),dv]~period+distance+period*distance + distance2+period*distance2, data=temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),])


	if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear squared"]<-NA
results[[j]]$se[results[[j]]$model.name=="loc. linear squared"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="loc. linear squared"]<-NA
		}
		
	if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear squared"]<-m$coefficients["period=post-memo"]
new<-robcov(m, cluster=temp[,clust.var])
results[[j]]$se[results[[j]]$model.name=="loc. linear squared"]<-sqrt(vcov(new)["period=post-memo","period=post-memo"] ) 
results[[j]]$se.std[results[[j]]$model.name=="loc. linear squared"]<-summary(m2)$coefficients["periodpost-memo","Std. Error"]
}

##local linear + lagged hit rate
m<-ols(temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),dv]~period+distance+period*distance + lag.hit.new, data=temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),],model=TRUE, x=TRUE, y=TRUE)

m2<-lm(temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),dv]~period+distance+period*distance + lag.hit.new, data=temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),])


	if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear + lagged dv"]<-NA
results[[j]]$se[results[[j]]$model.name=="loc. linear + lagged dv"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="loc. linear + lagged dv"]<-NA
}

	if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear + lagged dv"]<-m$coefficients["period=post-memo"]
new<-robcov(m, cluster=temp[,clust.var])
results[[j]]$se[results[[j]]$model.name=="loc. linear + lagged dv"]<-sqrt(vcov(new)["period=post-memo","period=post-memo"] ) 
results[[j]]$se.std[results[[j]]$model.name=="loc. linear + lagged dv"]<-summary(m2)$coefficients["periodpost-memo","Std. Error"]
}



##local linear squared + lagged hit rate
m<-ols(temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),dv]~period+distance+period*distance + distance2+period*distance2 + lag.hit.new, data=temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),],model=TRUE, x=TRUE, y=TRUE)

m2<-lm(temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),dv]~period+distance+period*distance + distance2+period*distance2 + lag.hit.new, data=temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),])


	if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="squared + lagged dv"]<-NA
results[[j]]$se[results[[j]]$model.name=="squared + lagged dv"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="squared + lagged dv"]<-NA}


	if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="squared + lagged dv"]<-m$coefficients["period=post-memo"]
new<-robcov(m, cluster=temp[,clust.var])
results[[j]]$se[results[[j]]$model.name=="squared + lagged dv"]<-sqrt(vcov(new)["period=post-memo","period=post-memo"] ) 
results[[j]]$se.std[results[[j]]$model.name=="squared + lagged dv"]<-summary(m2)$coefficients["periodpost-memo","Std. Error"]
}



##cubic
m<-ols(temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),dv]~period+distance+period*distance + distance2+period*distance2 +  distance3+period*distance3, data=temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),],model=TRUE, x=TRUE, y=TRUE)
m2<-lm(temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),dv]~period+distance+period*distance + distance2+period*distance2 +  distance3+period*distance3, data=temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),])


	if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="cubic"]<-NA
results[[j]]$se[results[[j]]$model.name=="cubic"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="cubic"]<-NA
}

	if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="cubic"]<-m$coefficients["period=post-memo"]
new<-robcov(m, cluster=temp[,clust.var])
results[[j]]$se[results[[j]]$model.name=="cubic"]<-sqrt(vcov(new)["period=post-memo","period=post-memo"] ) 
results[[j]]$se.std[results[[j]]$model.name=="cubic"]<-summary(m2)$coefficients["periodpost-memo","Std. Error"]
}

##difference in means
m<-ols(temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),dv]~period, data=temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),],model=TRUE, x=TRUE, y=TRUE)
m2<-lm(temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),dv]~period, data=temp[temp$dayno>=(cut2-N[j]) & temp$dayno<(cut2+N[j]),])

new<-robcov(m, cluster=temp[,clust.var])
results[[j]]$coef[results[[j]]$model.name=="diff"]<-m$coefficients["period=post-memo"]
results[[j]]$se[results[[j]]$model.name=="diff"]<-sqrt(vcov(new)["period=post-memo","period=post-memo"] ) 
results[[j]]$se.std[results[[j]]$model.name=="diff"]<-summary(m2)$coefficients["periodpost-memo","Std. Error"]

}
}

}

if(cluster!=T|is.null(cluster)){
for(j in 1:length(N)){
	
		data<-data[order(data$distance), ]

	
##local linear
m<-lm(data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),dv]~period+distance+period*distance, data=data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),])


if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear"]<-NA
results[[j]]$se[results[[j]]$model.name=="loc. linear"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="loc. linear"]<-NA
	}

if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear"]<-summary(m)$coefficients["periodpost-memo","Estimate"]
new<-coeftest(m, vcov.=vcovHAC(m, order.by=m$model$distance))
results[[j]]$se[results[[j]]$model.name=="loc. linear"]<-new["periodpost-memo","Std. Error"]
results[[j]]$se.std[results[[j]]$model.name=="loc. linear"]<-summary(m)$coefficients["periodpost-memo","Std. Error"]
}

##local linear squared
m<-lm(data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),dv]~period+distance+period*distance + I(distance^2)+period*I(distance^2), data=data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),])

	if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear squared"]<-NA
results[[j]]$se[results[[j]]$model.name=="loc. linear squared"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="loc. linear squared"]<-NA
		}
		
	if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear squared"]<-summary(m)$coefficients["periodpost-memo","Estimate"]
new<-coeftest(m, vcov.=vcovHAC(m, order.by=m$model$distance))
results[[j]]$se[results[[j]]$model.name=="loc. linear squared"]<-new["periodpost-memo","Std. Error"]
results[[j]]$se.std[results[[j]]$model.name=="loc. linear squared"]<-summary(m)$coefficients["periodpost-memo","Std. Error"]
}

##local linear + lagged hit rate
m<-lm(data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),dv]~period+distance+period*distance + lag.hit.new, data=data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),])


	if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear + lagged dv"]<-NA
results[[j]]$se[results[[j]]$model.name=="loc. linear + lagged dv"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="loc. linear + lagged dv"]<-NA
}

	if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear + lagged dv"]<-summary(m)$coefficients["periodpost-memo","Estimate"]
new<-coeftest(m, vcov.=vcovHAC(m, order.by=m$model$distance))
results[[j]]$se[results[[j]]$model.name=="loc. linear + lagged dv"]<-new["periodpost-memo","Std. Error"]
results[[j]]$se.std[results[[j]]$model.name=="loc. linear + lagged dv"]<-summary(m)$coefficients["periodpost-memo","Std. Error"]
}



##local linear squared + lagged hit rate
m<-lm(data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),dv]~period+distance+period*distance + I(distance^2)+period*I(distance^2) + lag.hit.new, data=data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),])


	if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="squared + lagged dv"]<-NA
results[[j]]$se[results[[j]]$model.name=="squared + lagged dv"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="squared + lagged dv"]<-NA}


	if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="squared + lagged dv"]<-summary(m)$coefficients["periodpost-memo","Estimate"]
new<-coeftest(m, vcov.=vcovHAC(m, order.by=m$model$distance))
results[[j]]$se[results[[j]]$model.name=="squared + lagged dv"]<-new["periodpost-memo","Std. Error"]
results[[j]]$se.std[results[[j]]$model.name=="squared + lagged dv"]<-summary(m)$coefficients["periodpost-memo","Std. Error"]}



##cubic
m<-lm(data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),dv]~period+distance+period*distance + I(distance^2)+period*I(distance^2) +  I(distance^3)+period*I(distance^3), data=data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),])

	if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="cubic"]<-NA
results[[j]]$se[results[[j]]$model.name=="cubic"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="cubic"]<-NA
}

	if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="cubic"]<-summary(m)$coefficients["periodpost-memo","Estimate"]
new<-coeftest(m, vcov.=vcovHAC(m, order.by=m$model$distance))
results[[j]]$se[results[[j]]$model.name=="cubic"]<-new["periodpost-memo","Std. Error"]
results[[j]]$se.std[results[[j]]$model.name=="cubic"]<-summary(m)$coefficients["periodpost-memo","Std. Error"]
}

##difference in means
temp<-data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),]
temp<-temp[order(temp$dayno),]
m<-lm(temp[,dv]~period, data=temp)
new<-coeftest(m, vcov.=vcovHAC(m))
results[[j]]$coef[results[[j]]$model.name=="diff"]<-summary(m)$coefficients["periodpost-memo","Estimate"]
results[[j]]$se[results[[j]]$model.name=="diff"]<-new["periodpost-memo","Std. Error"]
results[[j]]$se.std[results[[j]]$model.name=="diff"]<-summary(m)$coefficients["periodpost-memo","Std. Error"]

}
}


loc.linear.coef<-NA
loc.linear.squared.coef<-NA
cubic.coef<-NA
diff.coef<-NA
lagged.coef<-NA
lagged.squared.coef<-NA

##standard errors
loc.linear.squared.se<-NA
lagged.se<-NA
lagged.squared.se<-NA
loc.linear.se<-NA
diff.se<-NA
cubic.se<-NA

loc.linear.squared.se.std<-NA
lagged.se.std<-NA
lagged.squared.se.std<-NA
loc.linear.se.std<-NA
diff.se.std<-NA
cubic.se.std<-NA

for(i in min(N):length(results)){
	
	loc.linear.coef<-c(loc.linear.coef, results[[i]]$coef[results[[i]]$model.name=="loc. linear"])
	loc.linear.se<-c(loc.linear.se, results[[i]]$se[results[[i]]$model.name=="loc. linear"])
	loc.linear.se.std<-c(loc.linear.se.std, results[[i]]$se.std[results[[i]]$model.name=="loc. linear"])
	
	loc.linear.squared.coef<-c(loc.linear.squared.coef, results[[i]]$coef[results[[i]]$model.name=="loc. linear squared"])
	loc.linear.squared.se<-c(loc.linear.squared.se, results[[i]]$se[results[[i]]$model.name=="loc. linear squared"])
	loc.linear.squared.se.std<-c(loc.linear.squared.se.std, results[[i]]$se.std[results[[i]]$model.name=="loc. linear squared"])
	
	lagged.coef<-c(lagged.coef, results[[i]]$coef[results[[i]]$model.name=="loc. linear + lagged dv"])
	lagged.se<-c(lagged.se, results[[i]]$se[results[[i]]$model.name=="loc. linear + lagged dv"])
	lagged.se.std<-c(lagged.se.std, results[[i]]$se.std[results[[i]]$model.name=="loc. linear + lagged dv"])
	
	lagged.squared.coef<-c(lagged.squared.coef, results[[i]]$coef[results[[i]]$model.name=="squared + lagged dv"])
	lagged.squared.se<-c(lagged.squared.se, results[[i]]$se[results[[i]]$model.name=="squared + lagged dv"])
	lagged.squared.se.std<-c(lagged.squared.se.std, results[[i]]$se.std[results[[i]]$model.name=="squared + lagged dv"])
	
	
	cubic.coef<-c(cubic.coef, results[[i]]$coef[results[[i]]$model.name=="cubic"])
	cubic.se<-c(cubic.se, results[[i]]$se[results[[i]]$model.name=="cubic"])
	cubic.se.std<-c(cubic.se.std, results[[i]]$se.std[results[[i]]$model.name=="cubic"])

	diff.coef<-c(diff.coef, results[[i]]$coef[results[[i]]$model.name=="diff"])
	diff.se<-c(diff.se, results[[i]]$se[results[[i]]$model.name=="diff"])	
	diff.se.std<-c(diff.se.std, results[[i]]$se.std[results[[i]]$model.name=="diff"])	

}


print(length(loc.linear.coef))
print(length(loc.linear.se))


##organize results
if(cluster==T){
	
m1<-cbind.data.frame(coef=loc.linear.coef, se=loc.linear.se, se.std=loc.linear.se.std, name=rep("loc.linear",length(loc.linear.coef)))

m1$se.larger<-NA
for(i in 1:length(m1$se)){
	if( !is.na(m1$se[i]>m1$se.std[i]) ){
if( m1$se[i]>m1$se.std[i] ){m1$se.larger[i]<-m1$se[i]}
if( m1$se[i]<m1$se.std[i] ){m1$se.larger[i]<-m1$se.std[i]}}
}
m1<-m1[,c("coef","se.larger","name")]
colnames(m1)<-c("coef","se","name")
m1$lb<-m1$coef-1.96*m1$se
m1$ub<-m1$coef+1.96*m1$se


m2<-cbind.data.frame(coef= loc.linear.squared.coef, se=loc.linear.squared.se,se.std=loc.linear.squared.se.std,  name=rep("loc.linear.squared",length(loc.linear.squared.coef)))
m2$se.larger<-NA
for(i in 1:length(m2$se)){
	if( !is.na(m2$se[i]>m2$se.std[i]) ){
if(m2$se[i]>m2$se.std[i]){m2$se.larger[i]<-m2$se[i]}
if(m2$se[i]<m2$se.std[i]){m2$se.larger[i]<-m2$se.std[i]}}
}
m2<-m2[,c("coef","se.larger","name")]
colnames(m2)<-c("coef","se","name")
m2$lb<-m2$coef-1.96*m2$se
m2$ub<-m2$coef+1.96*m2$se

m3<-cbind.data.frame(coef= lagged.coef, se=lagged.se,se.std=lagged.se.std, name=rep("lagged",length(lagged.coef)))
m3$se.larger<-NA
for(i in 1:length(m3$se)){
	if( !is.na(m3$se[i]>m3$se.std[i]) ){
if(m3$se[i]>m3$se.std[i]){m3$se.larger[i]<-m3$se[i]}
if(m3$se[i]<m3$se.std[i]){m3$se.larger[i]<-m3$se.std[i]}}
}
m3<-m3[,c("coef","se.larger","name")]
colnames(m3)<-c("coef","se","name")
m3$lb<-m3$coef-1.96*m3$se
m3$ub<-m3$coef+1.96*m3$se


m4<-cbind.data.frame(coef= lagged.squared.coef, se=lagged.squared.se,se.std=lagged.squared.se.std, name=rep("lagged.squared",length(lagged.squared.coef)))
m4$se.larger<-NA
for(i in 1:length(m4$se)){
	if( !is.na(m4$se[i]>m4$se.std[i]) ){
if(m4$se[i]>m4$se.std[i]){m4$se.larger[i]<-m4$se[i]}
if(m4$se[i]<m4$se.std[i]){m4$se.larger[i]<-m4$se.std[i]}}
}
m4<-m4[,c("coef","se.larger","name")]
colnames(m4)<-c("coef","se","name")
m4$lb<-m4$coef-1.96*m4$se
m4$ub<-m4$coef+1.96*m4$se


m5<-cbind.data.frame(coef= diff.coef,se=diff.se,se.std=diff.se.std, name=rep("diff in means",length(diff.coef)))
m5$se.larger<-NA
for(i in 1:length(m5$se)){
	if( !is.na(m5$se[i]>m5$se.std[i]) ){
if(m5$se[i]>m5$se.std[i]){m5$se.larger[i]<-m5$se[i]}
if(m5$se[i]<m5$se.std[i]){m5$se.larger[i]<-m5$se.std[i]}}
}
m5<-m5[,c("coef","se.larger","name")]
colnames(m5)<-c("coef","se","name")
m5$lb<-m5$coef-1.96*m5$se
m5$ub<-m5$coef+1.96*m5$se



m6<-cbind.data.frame(coef= cubic.coef,se=cubic.se,se.std=cubic.se.std, name=rep("cubic",length(cubic.coef)))
m6$se.larger<-NA
for(i in 1:length(m6$se)){
	if( !is.na(m6$se[i]>m6$se.std[i]) ){
if(m6$se[i]>m6$se.std[i]){m6$se.larger[i]<-m6$se[i]}
if(m6$se[i]<m6$se.std[i]){m6$se.larger[i]<-m6$se.std[i]}}
}
m6<-m6[,c("coef","se.larger","name")]
colnames(m6)<-c("coef","se","name")
m6$lb<-m6$coef-1.96*m6$se
m6$ub<-m6$coef+1.96*m6$se

d<-list(results,m1,m2,m3,m4,m5,m6)
names(d)<-c("results","m1","m2","m3","m4","m5","m6")


d[["m1"]]<-d[["m1"]][-1,]
d[["m2"]]<-d[["m2"]][-1,]
d[["m3"]]<-d[["m3"]][-1,]
d[["m4"]]<-d[["m4"]][-1,]
d[["m5"]]<-d[["m5"]][-1,]
d[["m6"]]<-d[["m6"]][-1,]



	
}




#organize results
if(cluster!=T|is.null(cluster)){
m1<-cbind.data.frame(coef=loc.linear.coef, se=loc.linear.se, se.std=loc.linear.se.std, name=rep("loc.linear",max(N)+1))

m1$se.larger<-NA
for(i in 1:length(m1$se)){
	if( !is.na(m1$se[i]>m1$se.std[i]) ){
if( m1$se[i]>m1$se.std[i] ){m1$se.larger[i]<-m1$se[i]}
if( m1$se[i]<m1$se.std[i] ){m1$se.larger[i]<-m1$se.std[i]}}
}
m1<-m1[,c("coef","se.larger","name")]
colnames(m1)<-c("coef","se","name")
m1$lb<-m1$coef-1.96*m1$se
m1$ub<-m1$coef+1.96*m1$se


m2<-cbind.data.frame(coef= loc.linear.squared.coef, se=loc.linear.squared.se,se.std=loc.linear.squared.se.std,  name=rep("loc.linear.squared",max(N)+1))
m2$se.larger<-NA
for(i in 1:length(m2$se)){
	if( !is.na(m2$se[i]>m2$se.std[i]) ){
if(m2$se[i]>m2$se.std[i]){m2$se.larger[i]<-m2$se[i]}
if(m2$se[i]<m2$se.std[i]){m2$se.larger[i]<-m2$se.std[i]}}
}
m2<-m2[,c("coef","se.larger","name")]
colnames(m2)<-c("coef","se","name")
m2$lb<-m2$coef-1.96*m2$se
m2$ub<-m2$coef+1.96*m2$se

m3<-cbind.data.frame(coef= lagged.coef, se=lagged.se,se.std=lagged.se.std, name=rep("lagged",max(N)+1))
m3$se.larger<-NA
for(i in 1:length(m3$se)){
	if( !is.na(m3$se[i]>m3$se.std[i]) ){
if(m3$se[i]>m3$se.std[i]){m3$se.larger[i]<-m3$se[i]}
if(m3$se[i]<m3$se.std[i]){m3$se.larger[i]<-m3$se.std[i]}}
}
m3<-m3[,c("coef","se.larger","name")]
colnames(m3)<-c("coef","se","name")
m3$lb<-m3$coef-1.96*m3$se
m3$ub<-m3$coef+1.96*m3$se


m4<-cbind.data.frame(coef= lagged.squared.coef, se=lagged.squared.se,se.std=lagged.squared.se.std, name=rep("lagged.squared",max(N)+1))
m4$se.larger<-NA
for(i in 1:length(m4$se)){
	if( !is.na(m4$se[i]>m4$se.std[i]) ){
if(m4$se[i]>m4$se.std[i]){m4$se.larger[i]<-m4$se[i]}
if(m4$se[i]<m4$se.std[i]){m4$se.larger[i]<-m4$se.std[i]}}
}
m4<-m4[,c("coef","se.larger","name")]
colnames(m4)<-c("coef","se","name")
m4$lb<-m4$coef-1.96*m4$se
m4$ub<-m4$coef+1.96*m4$se


m5<-cbind.data.frame(coef= diff.coef,se=diff.se,se.std=diff.se.std, name=rep("diff in means",max(N)+1))
m5$se.larger<-NA
for(i in 1:length(m5$se)){
	if( !is.na(m5$se[i]>m5$se.std[i]) ){
if(m5$se[i]>m5$se.std[i]){m5$se.larger[i]<-m5$se[i]}
if(m5$se[i]<m5$se.std[i]){m5$se.larger[i]<-m5$se.std[i]}}
}
m5<-m5[,c("coef","se.larger","name")]
colnames(m5)<-c("coef","se","name")
m5$lb<-m5$coef-1.96*m5$se
m5$ub<-m5$coef+1.96*m5$se



m6<-cbind.data.frame(coef= cubic.coef,se=cubic.se,se.std=cubic.se.std, name=rep("cubic",max(N)+1))
m6$se.larger<-NA
for(i in 1:length(m6$se)){
	if( !is.na(m6$se[i]>m6$se.std[i]) ){
if(m6$se[i]>m6$se.std[i]){m6$se.larger[i]<-m6$se[i]}
if(m6$se[i]<m6$se.std[i]){m6$se.larger[i]<-m6$se.std[i]}}
}
m6<-m6[,c("coef","se.larger","name")]
colnames(m6)<-c("coef","se","name")
m6$lb<-m6$coef-1.96*m6$se
m6$ub<-m6$coef+1.96*m6$se
}


d<-list(results,m1,m2,m3,m4,m5,m6)
names(d)<-c("results","m1","m2","m3","m4","m5","m6")


d[["m1"]]<-d[["m1"]][-1,]
d[["m2"]]<-d[["m2"]][-1,]
d[["m3"]]<-d[["m3"]][-1,]
d[["m4"]]<-d[["m4"]][-1,]
d[["m5"]]<-d[["m5"]][-1,]
d[["m6"]]<-d[["m6"]][-1,]


return(d)
}







if(model=="logit"){
	
		cut2<-memo

	
	form<-as.formula(paste(dv,"~","dayno"))

means<-summaryBy(form, data=data, FUN=mean, na.rm=T)##get daily hit rates
means$lag.hit<-NA
means<-means[order(means$dayno, decreasing=F ),]




for(i in 2:length(means$dayno)){
	
	means$lag.hit[i]<-means[ i-1  ,   2]
	
}

names(means)[3]<-"lag.hit.new"


data<-merge(data, means, by="dayno")
head(data)
		


#make storage object
##make nested lists
diffs.list.dm<-list(NA)
diffs.list.l<-list(NA)
diffs.list.ll<-list(NA)
diffs.list.s<-list(NA)
diffs.list.sl<-list(NA)
diffs.list.c<-list(NA)


logit.boot<-list(NA)
lag<-data$lag.hit[data$dayno==(cut2)][1]

#start bootstrap loop
for(j in 1:length(N)){
	
	#get data set of bandwidth N
sw2<-data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),c("weapon.hit","weekday","period","distance","month","year","lag.hit","dayno")]
dim(sw2)


	
##get bootstrapped data set indices at bandwidth N
ind<-list(NA)##indices for bootstrapped data set at bandwidth N
rownames(sw2)<-1:nrow(sw2)
head(rownames(sw2))
set.seed(31415)

for(k in 1:1000){
	ind[[k]]<-sample(x=1:nrow(sw2),size=nrow(sw2),replace=T)##get bootstrapped sample indices
}	


diffs.list.dm[[j]]<-NA
##diff in means model
for(b in 1:1000){
	##make storage list for this model at bandwidth N
m<-glm(weapon.hit ~ period, data=sw2[ind[[b]],], family=binomial(link=logit))
	d2<-predict(m, newdata=cbind.data.frame(period="post-memo"), type="response")
	d1<-predict(m, newdata=cbind.data.frame(period="pre-memo"), type="response")
	diffs.list.dm[[j]][b]<-d2-d1
	
	#names(diffs.list.dm[[j]])<-paste("Diff in Means at bandwidth",N, sep=""  )
	
	if(b==1000){print(paste("Diff in Means at bandwidth ",N[j], sep=""  ))}
names(diffs.list.dm[[j]])<-paste("Diff in means bandwidth ",N[j], sep=""  )



}##close bootstrap loop at bandwidth N

logit.boot[[1]]<-diffs.list.dm
save(logit.boot, file="diffs.list.local.logit.Rdata")


diffs.list.l[[j]]<-NA 
##linear
for(b in 1:1000){
m<-glm(weapon.hit ~ period+distance+period*distance, data=sw2[ind[[b]],], family=binomial(link=logit))
		
	if(length(table(is.na(m$coefficients)))>1){
	diffs.list.ll[[j]][b]<-NA}
	
	if(length(table(is.na(m$coefficients)))==1){
	d2<-predict(m, newdata=cbind.data.frame(period="post-memo", distance=0), type="response")
	d1<-predict(m, newdata=cbind.data.frame(period="pre-memo",distance=0), type="response")
	diffs.list.l[[j]][b]<-d2-d1}
	
		if(b==1000){print(paste("Linear at bandwidth ",N[j], sep=""  ))}
		names(diffs.list.l[[j]])<-paste("Linear at bandwidth ",N[j], sep=""  )

	

}##close bootstrap loop at bandwidth N


logit.boot[[2]]<-diffs.list.l
save(logit.boot, file="diffs.list.local.logit.Rdata")



##linear + lagged DV
diffs.list.ll[[j]]<-NA
for(b in 1:1000){
m<-glm(weapon.hit ~ period+distance+period*distance+lag.hit, data=sw2[ind[[b]],], family=binomial(link=logit))
	d2<-predict(m, newdata=cbind.data.frame(period="post-memo", distance=0, lag.hit=lag), type="response")
	d1<-predict(m, newdata=cbind.data.frame(period="pre-memo", distance=0, lag.hit=lag), type="response")
	
	if(length(table(is.na(m$coefficients)))>1){
	diffs.list.ll[[j]][b]<-NA}
	if(length(table(is.na(m$coefficients)))==1){
	diffs.list.ll[[j]][b]<-d2-d1}	

	if(b==1000){print(paste("Linear + controls at bandwidth ",N[j], sep=""  ))}
names(diffs.list.ll[[j]])<-paste("Linear + controls at bandwidth ",N[j], sep=""  )


}##close bootstrap loop at bandwidth N

logit.boot[[3]]<-diffs.list.ll
save(logit.boot, file="diffs.list.local.logit.Rdata")




##squared
diffs.list.s[[j]]<-NA

for(b in 1:1000){
m<-glm(weapon.hit ~ period+distance+period*distance+I(distance^2)+period*I(distance^2), data=sw2[ind[[b]],], family=binomial(link=logit))

	if(length(table(is.na(m$coefficients)))>1){
	diffs.list.s[[j]][b]<-NA}
	

	if(length(table(is.na(m$coefficients)))==1){
			d2<-predict(m, newdata=cbind.data.frame(period="post-memo",distance=0), type="response")
	d1<-predict(m, newdata=cbind.data.frame(period="pre-memo", distance=0), type="response")
	diffs.list.s[[j]][b]<-d2-d1}		

	if(b==1000){print(paste("Squared at bandwidth ",N[j], sep=""  ))}
names(diffs.list.s[[j]])<-paste("Squared at bandwidth ",N[j], sep=""  )


}##close bootstrap loop at bandwidth N

logit.boot[[4]]<-diffs.list.s
save(logit.boot, file="diffs.list.local.logit.Rdata")




##squared + lagged dv
diffs.list.sl[[j]]<-NA 
for(b in 1:1000){
m<-glm(weapon.hit ~ period+distance+period*distance+I(distance^2)+period*I(distance^2)+lag.hit, data=sw2[ind[[b]],], family=binomial(link=logit))

	if(length(table(is.na(m$coefficients)))>1){
	diffs.list.sl[[j]][b]<-NA}
	

	

	if(length(table(is.na(m$coefficients)))==1){
			d2<-predict(m, newdata=cbind.data.frame(period="post-memo",distance=0, lag.hit=lag), type="response")
	d1<-predict(m, newdata=cbind.data.frame(period="pre-memo", distance=0, lag.hit=lag), type="response")
	diffs.list.sl[[j]][b]<-d2-d1}		

	

	if(b==1000){print(paste("Squared + controls at bandwidth ",N[j], sep=""  ))}
names(diffs.list.sl[[j]])<-paste("Squared + controls at bandwidth ",N[j], sep=""  )


}##close bootstrap loop at bandwidth N

logit.boot[[5]]<-diffs.list.sl
save(logit.boot, file="diffs.list.local.logit.Rdata")



##cubic
diffs.list.c[[j]]<-NA
for(b in 1:1000){
m<-glm(weapon.hit ~ period+distance+period*distance+I(distance^2)+period*I(distance^2)+I(distance^3)+period*I(distance^3), data=sw2[ind[[b]],], family=binomial(link=logit))

	if(length(table(is.na(m$coefficients)))>1){
	diffs.list.c[[j]][b]<-NA}


	if(length(table(is.na(m$coefficients)))==1){
	d2<-predict(m,  newdata=cbind.data.frame(period="post-memo",distance=0), type="response")
	d1<-predict(m,  newdata=cbind.data.frame(period="pre-memo",distance=0), type="response")
	
	diffs.list.c[[j]][b]<-d2-d1}		

if(b==1000)	{print(paste("cubic at bandwidth ",N[j], sep=""  ))}
names(diffs.list.c[[j]])<-paste("cubic at bandwidth ",N[j], sep=""  )


}##close bootstrap loop at bandwidth N

logit.boot[[6]]<-diffs.list.c
save(logit.boot, file="diffs.list.local.logit.Rdata")


}#close bandwidth loop

return(logit.boot)
}#close logit loop







##interaction models
if(model%in%c("ols",NULL)  & int=="yes"){
	
	if(placebo==T){
	cut2<-memo
	
data$period<-NA
data$period[data$dayno<cut2]<-"pre-memo"
data$period[data$dayno>=cut2 ]<-"post-memo"
data$period<-factor(data$period, levels=c("pre-memo", "post-memo"))
table(data$period)
	
	data$distance<-data$dayno-cut2
	table(data$distance)

	}
	
form<-as.formula(paste(dv,"~","dayno"))

means<-summaryBy(form, data=data, FUN=mean, na.rm=T)##get daily hit rates
means$lag.hit<-NA
means<-means[order(means$dayno, decreasing=F ),]




for(i in 2:length(means$dayno)){
	
	means$lag.hit[i]<-means[ i-1  ,   2]
	
}

names(means)[3]<-"lag.hit.new"


data<-merge(data, means, by="dayno")
head(data)
		
		data<-data[order(data$distance), ]

	
model.names<-c("loc. linear","loc. linear squared","loc. linear + lagged dv", "squared + lagged dv","cubic","diff")
results<-list()

coefs<-NA
ses<-NA
se.std<-NA


for(i in 1:length(N)){
	
	results[[i]]<-cbind.data.frame(model.name=model.names, coef=coefs, se=ses, se.std=se.std)
}


names(results)<-N


data$int<-data[,int.var]


for(j in 1:length(N)){
	
	
	
##local linear
m<-lm(data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),dv]~period+distance+period*distance+int+period*int+distance*int+period*int*distance, data=data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),])


if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear"]<-NA
results[[j]]$se[results[[j]]$model.name=="loc. linear"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="loc. linear"]<-NA
	}

if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear"]<-summary(m)$coefficients["periodpost-memo:int","Estimate"]
new<-coeftest(m, vcov.=vcovHAC(m, order.by=m$model$distance))
results[[j]]$se[results[[j]]$model.name=="loc. linear"]<-new["periodpost-memo:int","Std. Error"]
results[[j]]$se.std[results[[j]]$model.name=="loc. linear"]<-summary(m)$coefficients["periodpost-memo:int","Std. Error"]
}

##local linear squared
m<-lm(data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),dv]~period+distance+period*distance + I(distance^2)+period*I(distance^2)+int+period*int+distance*int+I(distance^2)*int+period*int*distance+period*int*I(distance^2), data=data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),])

	if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear squared"]<-NA
results[[j]]$se[results[[j]]$model.name=="loc. linear squared"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="loc. linear squared"]<-NA
		}
		
	if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear squared"]<-summary(m)$coefficients["periodpost-memo:int","Estimate"]
new<-coeftest(m, vcov.=vcovHAC(m, order.by=m$model$distance))
results[[j]]$se[results[[j]]$model.name=="loc. linear squared"]<-new["periodpost-memo:int","Std. Error"]
results[[j]]$se.std[results[[j]]$model.name=="loc. linear squared"]<-summary(m)$coefficients["periodpost-memo:int","Std. Error"]
}

##local linear + lagged hit rate
m<-lm(data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),dv]~period+distance+period*distance + lag.hit.new+int+period*int+distance*int+period*distance*int, data=data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),])


	if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear + lagged dv"]<-NA
results[[j]]$se[results[[j]]$model.name=="loc. linear + lagged dv"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="loc. linear + lagged dv"]<-NA
}

	if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="loc. linear + lagged dv"]<-summary(m)$coefficients["periodpost-memo:int","Estimate"]
new<-coeftest(m, vcov.=vcovHAC(m, order.by=m$model$distance))
results[[j]]$se[results[[j]]$model.name=="loc. linear + lagged dv"]<-new["periodpost-memo:int","Std. Error"]
results[[j]]$se.std[results[[j]]$model.name=="loc. linear + lagged dv"]<-summary(m)$coefficients["periodpost-memo:int","Std. Error"]
}



##local linear squared + lagged hit rate
m<-lm(data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),dv]~period+distance+period*distance + I(distance^2)+period*I(distance^2) + lag.hit.new +int+period*int+distance*int+I(distance^2)*int+period*int*distance+period*int*I(distance^2), data=data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),])


	if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="squared + lagged dv"]<-NA
results[[j]]$se[results[[j]]$model.name=="squared + lagged dv"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="squared + lagged dv"]<-NA}


	if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="squared + lagged dv"]<-summary(m)$coefficients["periodpost-memo:int","Estimate"]
new<-coeftest(m, vcov.=vcovHAC(m), order.by=m$model$distance)
results[[j]]$se[results[[j]]$model.name=="squared + lagged dv"]<-new["periodpost-memo:int","Std. Error"]
results[[j]]$se.std[results[[j]]$model.name=="squared + lagged dv"]<-summary(m)$coefficients["periodpost-memo:int","Std. Error"]}



##cubic
m<-lm(data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),dv]~period+distance+period*distance + I(distance^2)+period*I(distance^2) +  I(distance^3)+period*I(distance^3)+int+period*int+distance*int+I(distance^2)*int+period*int*distance+period*int*I(distance^2)+int*I(distance^3)+int*period*I(distance^3), data=data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),])

	if(length(table(is.na(m$coefficients)))>1){
results[[j]]$coef[results[[j]]$model.name=="cubic"]<-NA
results[[j]]$se[results[[j]]$model.name=="cubic"]<-NA
results[[j]]$se.std[results[[j]]$model.name=="cubic"]<-NA
}

	if(length(table(is.na(m$coefficients)))==1){
results[[j]]$coef[results[[j]]$model.name=="cubic"]<-summary(m)$coefficients["periodpost-memo:int","Estimate"]
new<-coeftest(m, vcov.=vcovHAC(m, order.by=m$model$distance))
results[[j]]$se[results[[j]]$model.name=="cubic"]<-new["periodpost-memo:int","Std. Error"]
results[[j]]$se.std[results[[j]]$model.name=="cubic"]<-summary(m)$coefficients["periodpost-memo:int","Std. Error"]
}

##difference in means
temp<-data[data$dayno>=(cut2-N[j]) & data$dayno<(cut2+N[j]),]
temp<-temp[order(temp$distance),]

m<-lm(temp[,dv]~period+int+period*int, data=temp)
new<-coeftest(m, vcov.=vcovHAC(m))
results[[j]]$coef[results[[j]]$model.name=="diff"]<-summary(m)$coefficients["periodpost-memo:int","Estimate"]
results[[j]]$se[results[[j]]$model.name=="diff"]<-new["periodpost-memo:int","Std. Error"]
results[[j]]$se.std[results[[j]]$model.name=="diff"]<-summary(m)$coefficients["periodpost-memo:int","Std. Error"]

}

loc.linear.coef<-NA
loc.linear.squared.coef<-NA
cubic.coef<-NA
diff.coef<-NA
lagged.coef<-NA
lagged.squared.coef<-NA

##standard errors
loc.linear.squared.se<-NA
lagged.se<-NA
lagged.squared.se<-NA
loc.linear.se<-NA
diff.se<-NA
cubic.se<-NA

loc.linear.squared.se.std<-NA
lagged.se.std<-NA
lagged.squared.se.std<-NA
loc.linear.se.std<-NA
diff.se.std<-NA
cubic.se.std<-NA

for(i in 1:length(results)){
	
	loc.linear.coef<-c(loc.linear.coef, results[[i]]$coef[results[[i]]$model.name=="loc. linear"])
	loc.linear.se<-c(loc.linear.se, results[[i]]$se[results[[i]]$model.name=="loc. linear"])
	loc.linear.se.std<-c(loc.linear.se.std, results[[i]]$se.std[results[[i]]$model.name=="loc. linear"])
	
	loc.linear.squared.coef<-c(loc.linear.squared.coef, results[[i]]$coef[results[[i]]$model.name=="loc. linear squared"])
	loc.linear.squared.se<-c(loc.linear.squared.se, results[[i]]$se[results[[i]]$model.name=="loc. linear squared"])
	loc.linear.squared.se.std<-c(loc.linear.squared.se.std, results[[i]]$se.std[results[[i]]$model.name=="loc. linear squared"])
	
	lagged.coef<-c(lagged.coef, results[[i]]$coef[results[[i]]$model.name=="loc. linear + lagged dv"])
	lagged.se<-c(lagged.se, results[[i]]$se[results[[i]]$model.name=="loc. linear + lagged dv"])
	lagged.se.std<-c(lagged.se.std, results[[i]]$se.std[results[[i]]$model.name=="loc. linear + lagged dv"])
	
	lagged.squared.coef<-c(lagged.squared.coef, results[[i]]$coef[results[[i]]$model.name=="squared + lagged dv"])
	lagged.squared.se<-c(lagged.squared.se, results[[i]]$se[results[[i]]$model.name=="squared + lagged dv"])
	lagged.squared.se.std<-c(lagged.squared.se.std, results[[i]]$se.std[results[[i]]$model.name=="squared + lagged dv"])
	
	
	cubic.coef<-c(cubic.coef, results[[i]]$coef[results[[i]]$model.name=="cubic"])
	cubic.se<-c(cubic.se, results[[i]]$se[results[[i]]$model.name=="cubic"])
	cubic.se.std<-c(cubic.se.std, results[[i]]$se.std[results[[i]]$model.name=="cubic"])

	diff.coef<-c(diff.coef, results[[i]]$coef[results[[i]]$model.name=="diff"])
	diff.se<-c(diff.se, results[[i]]$se[results[[i]]$model.name=="diff"])	
	diff.se.std<-c(diff.se.std, results[[i]]$se.std[results[[i]]$model.name=="diff"])	

}


print(length(loc.linear.coef))
print(length(loc.linear.se))

m1<-cbind.data.frame(coef=loc.linear.coef, se=loc.linear.se, se.std=loc.linear.se.std, name=rep("loc.linear",max(N)+1))

m1$se.larger<-NA
for(i in 1:length(m1$se)){
	if( !is.na(m1$se[i]>m1$se.std[i]) ){
if( m1$se[i]>m1$se.std[i] ){m1$se.larger[i]<-m1$se[i]}
if( m1$se[i]<m1$se.std[i] ){m1$se.larger[i]<-m1$se.std[i]}}
}
m1<-m1[,c("coef","se.larger","name")]
colnames(m1)<-c("coef","se","name")
m1$lb<-m1$coef-1.96*m1$se
m1$ub<-m1$coef+1.96*m1$se


m2<-cbind.data.frame(coef= loc.linear.squared.coef, se=loc.linear.squared.se,se.std=loc.linear.squared.se.std,  name=rep("loc.linear.squared",max(N)+1))
m2$se.larger<-NA
for(i in 1:length(m2$se)){
	if( !is.na(m2$se[i]>m2$se.std[i]) ){
if(m2$se[i]>m2$se.std[i]){m2$se.larger[i]<-m2$se[i]}
if(m2$se[i]<m2$se.std[i]){m2$se.larger[i]<-m2$se.std[i]}}
}
m2<-m2[,c("coef","se.larger","name")]
colnames(m2)<-c("coef","se","name")
m2$lb<-m2$coef-1.96*m2$se
m2$ub<-m2$coef+1.96*m2$se

m3<-cbind.data.frame(coef= lagged.coef, se=lagged.se,se.std=lagged.se.std, name=rep("lagged",max(N)+1))
m3$se.larger<-NA
for(i in 1:length(m3$se)){
	if( !is.na(m3$se[i]>m3$se.std[i]) ){
if(m3$se[i]>m3$se.std[i]){m3$se.larger[i]<-m3$se[i]}
if(m3$se[i]<m3$se.std[i]){m3$se.larger[i]<-m3$se.std[i]}}
}
m3<-m3[,c("coef","se.larger","name")]
colnames(m3)<-c("coef","se","name")
m3$lb<-m3$coef-1.96*m3$se
m3$ub<-m3$coef+1.96*m3$se


m4<-cbind.data.frame(coef= lagged.squared.coef, se=lagged.squared.se,se.std=lagged.squared.se.std, name=rep("lagged.squared",max(N)+1))
m4$se.larger<-NA
for(i in 1:length(m4$se)){
	if( !is.na(m4$se[i]>m4$se.std[i]) ){
if(m4$se[i]>m4$se.std[i]){m4$se.larger[i]<-m4$se[i]}
if(m4$se[i]<m4$se.std[i]){m4$se.larger[i]<-m4$se.std[i]}}
}
m4<-m4[,c("coef","se.larger","name")]
colnames(m4)<-c("coef","se","name")
m4$lb<-m4$coef-1.96*m4$se
m4$ub<-m4$coef+1.96*m4$se


m5<-cbind.data.frame(coef= diff.coef,se=diff.se,se.std=diff.se.std, name=rep("diff in means",max(N)+1))
m5$se.larger<-NA
for(i in 1:length(m5$se)){
	if( !is.na(m5$se[i]>m5$se.std[i]) ){
if(m5$se[i]>m5$se.std[i]){m5$se.larger[i]<-m5$se[i]}
if(m5$se[i]<m5$se.std[i]){m5$se.larger[i]<-m5$se.std[i]}}
}
m5<-m5[,c("coef","se.larger","name")]
colnames(m5)<-c("coef","se","name")
m5$lb<-m5$coef-1.96*m5$se
m5$ub<-m5$coef+1.96*m5$se



m6<-cbind.data.frame(coef= cubic.coef,se=cubic.se,se.std=cubic.se.std, name=rep("cubic",max(N)+1))
m6$se.larger<-NA
for(i in 1:length(m6$se)){
	if( !is.na(m6$se[i]>m6$se.std[i]) ){
if(m6$se[i]>m6$se.std[i]){m6$se.larger[i]<-m6$se[i]}
if(m6$se[i]<m6$se.std[i]){m6$se.larger[i]<-m6$se.std[i]}}
}
m6<-m6[,c("coef","se.larger","name")]
colnames(m6)<-c("coef","se","name")
m6$lb<-m6$coef-1.96*m6$se
m6$ub<-m6$coef+1.96*m6$se



d<-list(results,m1,m2,m3,m4,m5,m6)
names(d)<-c("results","m1","m2","m3","m4","m5","m6")

#remove extra row of NA
d[["m1"]]<-d[["m1"]][-1,]
d[["m2"]]<-d[["m2"]][-1,]
d[["m3"]]<-d[["m3"]][-1,]
d[["m4"]]<-d[["m4"]][-1,]
d[["m5"]]<-d[["m5"]][-1,]
d[["m6"]]<-d[["m6"]][-1,]



return(d)
}








}#close function